/*
 #
 #  Files       : MatrixBase.h
 #                ( C++ header file )
 #
 #  Description : The SmallMatrix Library
 #                ( http://code.google.com/p/smallmatrix )
 #
 #  Copyright   : Olivier Juan
 #                ( http://www.mas.ecp.fr/vision/Personnel/juan/ )
 #
 #  License     : CeCILL-C
 #                ( http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html )
 #
 #  This software is governed by the CeCILL-C license under French law and
 #  abiding by the rules of distribution of free software.  You can  use,
 #  modify and or redistribute the software under the terms of the CeCILL-C
 #  license as circulated by CEA, CNRS and INRIA at the following URL
 #  "http://www.cecill.info".
 #
 #  As a counterpart to the access to the source code and  rights to copy,
 #  modify and redistribute granted by the license, users are provided only
 #  with a limited warranty  and the software's author,  the holder of the
 #  economic rights,  and the successive licensors  have only  limited
 #  liability.
 #
 #  In this respect, the user's attention is drawn to the risks associated
 #  with loading,  using,  modifying and/or developing or reproducing the
 #  software by the user in light of its specific status of free software,
 #  that may mean  that it is complicated to manipulate,  and  that  also
 #  therefore means  that it is reserved for developers  and  experienced
 #  professionals having in-depth computer knowledge. Users are therefore
 #  encouraged to load and test the software's suitability as regards their
 #  requirements in conditions enabling the security of their systems and/or
 #  data to be ensured and,  more generally, to use and operate it in the
 #  same conditions as regards security.
 #
 #  The fact that you are presently reading this means that you have had
 #  knowledge of the CeCILL-C license and that you accept its terms.
 #
*/

#ifndef _MATRIX_BASE_H
#define _MATRIX_BASE_H

#include <SmallMatrix/SmallMatrix.h>

NAMESPACE_BEGIN(SmallMatrix)

template <typename T,int SIZE,int M,int N> class Matrix_Base
{
public:
	typedef T self;
protected:
	reel Tab[SIZE];
	typename Matrix_Base<T,SIZE,M,N>::self& Child() {return static_cast<T&>(*this);}
	const typename Matrix_Base<T,SIZE,M,N>::self& Child() const {return static_cast<const T&>(*this);}
public:

	Matrix_Base() {
		Matrix_Base<T,SIZE,M,N>& ref=*this;
		for (int i=0;i<SIZE;i++)
			ref[i]=0;
	}
	Matrix_Base(const reel &v) {
		Matrix_Base<T,SIZE,M,N>& ref=*this;
		for (int i=0;i<SIZE;i++)
			ref[i]=v;
	}
	Matrix_Base(const reel tab[SIZE]) {
		Matrix_Base<T,SIZE,M,N>& ref=*this;
		for (int i=0;i<SIZE;i++)
			ref[i]=tab[i];
	}
	template <typename TT,int SSIZE> Matrix_Base(const Matrix_Base<TT,SSIZE,M,N>& copy) {
		typename Matrix_Strictly_Additive_Scheme_Trait<typename Matrix_Base<T,SIZE,M,N>::self ,typename Matrix_Base<TT,SSIZE,M,N>::self >::result_type& ref=Child(); 
		for (int i=0;i<M;i++)
			for (int j=0;j<N;j++)
				ref(i,j)=copy(i,j);
	}
	template <typename TT,int SSIZE> const typename Matrix_Base<TT,SSIZE,M,N>::self& operator=(const Matrix_Base<TT,SSIZE,M,N>& copy) {
		typename Matrix_Strictly_Additive_Scheme_Trait<typename Matrix_Base<T,SIZE,M,N>::self ,typename Matrix_Base<TT,SSIZE,M,N>::self >::result_type& ref=Child(); 
		for (int i=0;i<M;i++)
			for (int j=0;j<N;j++)
				ref(i,j)=copy(i,j);
		return copy;
	}

	Vecteur<M> col(int  j);
	Vecteur<N> row(int  i);
	reel _min();
	reel _min(int & idxmin);
	reel _max();
	reel _max(int & idxmax);

	inline typename Matrix_Base<T,SIZE,M,N>::self& sqrt();
	inline typename Matrix_Base<T,SIZE,M,N>::self operator-()const;
	inline typename Matrix_Base<T,SIZE,M,N>::self operator*(const reel& alpha)const;
	inline typename Matrix_Base<T,SIZE,M,N>::self operator/(const reel& alpha)const;
	inline typename Matrix_Base<T,SIZE,M,N>::self& operator*=(const reel& alpha);
	inline typename Matrix_Base<T,SIZE,M,N>::self& operator/=(const reel& alpha);
	inline typename Matrix_Base<T,SIZE,M,N>::self operator+(const reel& alpha)const;
	inline typename Matrix_Base<T,SIZE,M,N>::self operator-(const reel& alpha)const;
	inline typename Matrix_Base<T,SIZE,M,N>::self& operator+=(const reel& alpha);
	inline typename Matrix_Base<T,SIZE,M,N>::self& operator-=(const reel& alpha);
	inline reel& operator[](const int& idx);
	inline const reel& operator[](const int& idx)const;
	inline reel& operator()(const int& idxm,const int& idxn);
	inline const reel& operator()(const int& idxm,const int& idxn)const;
	inline typename Matrix_Base<T,SIZE,M,N>::self& set(const reel& value);
	inline reel norm()const;
	inline reel norm2()const;
	inline reel normL1()const;
	inline reel normL2()const;
	template <typename TT,int SSIZE> inline reel dot(const Matrix_Base<TT,SSIZE,M,N>& Two)const;
	template <typename TT,int SSIZE,int NN> inline typename Matrix_Product_Scheme_Trait< typename Matrix_Base<T,SIZE,M,N>::self, typename Matrix_Base<TT,SSIZE,N,NN>::self >::result_type operator*(const Matrix_Base<TT,SSIZE,N,NN>& M1) const;
	template <typename TT,int SSIZE> inline typename Matrix_Base<T,SIZE,M,N>::self& mult(const Matrix_Base<TT,SSIZE,M,N>& Two);
	template <typename TT,int SSIZE> inline typename Matrix_Base<T,SIZE,M,N>::self& div(const Matrix_Base<TT,SSIZE,M,N>& Two);
	template <typename TT,int SSIZE> inline bool operator==(const Matrix_Base<TT,SSIZE,M,N>& Two) const;
	template <typename TT,int SSIZE> inline typename Matrix_Additive_Scheme_Trait< typename Matrix_Base<T,SIZE,M,N>::self, typename Matrix_Base<TT,SSIZE,M,N>::self >::result_type operator+(const Matrix_Base<TT,SSIZE,M,N>& Two) const;
	template <typename TT,int SSIZE> inline typename Matrix_Additive_Scheme_Trait< typename Matrix_Base<T,SIZE,M,N>::self, typename Matrix_Base<TT,SSIZE,M,N>::self >::result_type operator-(const Matrix_Base<TT,SSIZE,M,N>& Two) const;
	template <typename TT,int SSIZE> inline typename Matrix_Base<T,SIZE,M,N>::self& operator+=(const Matrix_Base<TT,SSIZE,M,N>& Two);
	template <typename TT,int SSIZE> inline typename Matrix_Base<T,SIZE,M,N>::self& operator-=(const Matrix_Base<TT,SSIZE,M,N>& Two);
	inline void write(std::ostream& f) const;
	inline void read(std::istream& f);
};

template <typename T,int SIZE,int M,int N> Vecteur<M> Matrix_Base<T,SIZE,M,N>::col(int  j){
	typename Matrix_Base<T,SIZE,M,N>::self& ref=Child();
	Vecteur<M> _col;
	for(int i = 0 ; i < M ;i++)
		_col(i) = ref(i,j);
	return(_col);
}
template <typename T,int SIZE,int M,int N> Vecteur<N> Matrix_Base<T,SIZE,M,N>::row(int  i){
	typename Matrix_Base<T,SIZE,M,N>::self& ref=Child();
	Vecteur<N> _row;
	for(int j = 0 ; j < N ; j++)
		_row(j) = ref(i,j);
	return(_row);
}
template <typename T,int SIZE,int M,int N> reel Matrix_Base<T,SIZE,M,N>::_min(){
	reel sol = Tab[0];
	for(int i = 0 ; i < SIZE ; i++)
		if(Tab[i] < sol) sol = Tab[i];
	return(sol);
}
template <typename T,int SIZE,int M,int N> reel Matrix_Base<T,SIZE,M,N>::_min(int & idxmin){
	reel sol = Tab[0]; idxmin = 0;
	for(int i = 0 ; i < SIZE ; i++)
		if(Tab[i] < sol){
			sol = Tab[i];
			idxmin = i;
		}
	return(sol);
}
template <typename T,int SIZE,int M,int N> reel Matrix_Base<T,SIZE,M,N>::_max(){
	reel sol = Tab[0];
	for(int i = 0 ; i < SIZE ; i++)
		if(Tab[i] > sol) sol = Tab[i];
	return(sol);
}
template <typename T,int SIZE,int M,int N> reel Matrix_Base<T,SIZE,M,N>::_max(int & idxmax){
	reel sol = Tab[0]; idxmax = 0;
	for(int i = 0 ; i < SIZE ; i++)
		if(Tab[i] > sol){
			sol = Tab[i];
			idxmax = i;
		}
	return(sol);
}


template <typename T,int SIZE,int M,int N> inline typename Matrix_Base<T,SIZE,M,N>::self& Matrix_Base<T,SIZE,M,N>::sqrt()
{
	typename Matrix_Base<T,SIZE,M,N>::self& ref=Child();
	for (int i=0;i<SIZE;i++) {
		assert(ref[i]>=0);
		ref[i]=::sqrt(ref[i]);
	}
	return ref;
}
template <typename T,int SIZE,int M,int N> inline typename Matrix_Base<T,SIZE,M,N>::self Matrix_Base<T,SIZE,M,N>::operator-()const
{
	typename Matrix_Base<T,SIZE,M,N>::self temp(Child());
	for (int i=0;i<SIZE;i++)
		temp[i]*=-1;
	return temp;
}
template <typename T,int SIZE,int M,int N> inline typename Matrix_Base<T,SIZE,M,N>::self Matrix_Base<T,SIZE,M,N>::operator*(const reel& alpha)const
{
	typename Matrix_Base<T,SIZE,M,N>::self temp(Child());
	for (int i=0;i<SIZE;i++)
		temp[i]*=alpha;
	return temp;
}
template <typename T,int SIZE,int M,int N> inline typename Matrix_Base<T,SIZE,M,N>::self Matrix_Base<T,SIZE,M,N>::operator/(const reel& alpha)const
{
	assert(alpha!=0);
	typename Matrix_Base<T,SIZE,M,N>::self temp(Child());
	for (int i=0;i<SIZE;i++)
		temp[i]/=alpha;
	return temp;
}
template <typename T,int SIZE,int M,int N> inline typename Matrix_Base<T,SIZE,M,N>::self& Matrix_Base<T,SIZE,M,N>::operator*=(const reel& alpha)
{
	typename Matrix_Base<T,SIZE,M,N>::self& ref=Child();
	for (int i=0;i<SIZE;i++)
		ref[i]*=alpha;
	return ref;
}
template <typename T,int SIZE,int M,int N> inline typename Matrix_Base<T,SIZE,M,N>::self& Matrix_Base<T,SIZE,M,N>::operator/=(const reel& alpha)
{
	assert(alpha!=0);
	typename Matrix_Base<T,SIZE,M,N>::self& ref=Child();
	for (int i=0;i<SIZE;i++)
		ref[i]/=alpha;
	return ref;
}
template <typename T,int SIZE,int M,int N> inline typename Matrix_Base<T,SIZE,M,N>::self Matrix_Base<T,SIZE,M,N>::operator+(const reel& alpha)const
{
	typename Matrix_Base<T,SIZE,M,N>::self temp(Child());
	for (int i=0;i<SIZE;i++)
		temp[i]+=alpha;
	return temp;
}
template <typename T,int SIZE,int M,int N> inline typename Matrix_Base<T,SIZE,M,N>::self Matrix_Base<T,SIZE,M,N>::operator-(const reel& alpha)const
{
	typename Matrix_Base<T,SIZE,M,N>::self temp(Child());
	for (int i=0;i<SIZE;i++)
		temp[i]-=alpha;
	return temp;
}
template <typename T,int SIZE,int M,int N> inline typename Matrix_Base<T,SIZE,M,N>::self& Matrix_Base<T,SIZE,M,N>::operator+=(const reel& alpha)
{
	typename Matrix_Base<T,SIZE,M,N>::self& ref=Child();
	for (int i=0;i<SIZE;i++)
		ref[i]+=alpha;
	return ref;
}
template <typename T,int SIZE,int M,int N> inline typename Matrix_Base<T,SIZE,M,N>::self& Matrix_Base<T,SIZE,M,N>::operator-=(const reel& alpha)
{
	typename Matrix_Base<T,SIZE,M,N>::self& ref=Child();
	for (int i=0;i<SIZE;i++)
		ref[i]-=alpha;
	return ref;
}
template <typename T,int SIZE,int M,int N> inline reel& Matrix_Base<T,SIZE,M,N>::operator[](const int& idx)
{
	assert((idx>=0)&&(idx<SIZE));
	return Tab[idx];
}
template <typename T,int SIZE,int M,int N> inline const reel& Matrix_Base<T,SIZE,M,N>::operator[](const int& idx)const
{
	assert((idx>=0)&&(idx<SIZE));
	return Tab[idx];
}
template <typename T,int SIZE,int M,int N> inline reel& Matrix_Base<T,SIZE,M,N>::operator()(const int& idxm,const int& idxn)
{
	assert((idxm>=0)&&(idxm<M)&&(idxn>=0)&&(idxn<N));
	return Child()(idxm,idxn);
}
template <typename T,int SIZE,int M,int N> inline const reel& Matrix_Base<T,SIZE,M,N>::operator()(const int& idxm,const int& idxn)const
{
	assert((idxm>=0)&&(idxm<M)&&(idxn>=0)&&(idxn<N));
	return Child()(idxm,idxn);
}
template <typename T,int SIZE,int M,int N> inline typename Matrix_Base<T,SIZE,M,N>::self& Matrix_Base<T,SIZE,M,N>::set(const reel& value)
{
	typename Matrix_Base<T,SIZE,M,N>::self& ref=Child();
	for (int i=0;i<SIZE;i++)
		ref[i]=value;
	return ref;
}
template <typename T,int SIZE,int M,int N> inline reel Matrix_Base<T,SIZE,M,N>::norm()const
{
	return ::sqrt(norm2());
}
template <typename T,int SIZE,int M,int N> inline reel Matrix_Base<T,SIZE,M,N>::norm2()const
{
	return dot(*this);
}
template <typename T,int SIZE,int M,int N> inline reel Matrix_Base<T,SIZE,M,N>::normL1()const
{
	const Matrix_Base<T,SIZE,M,N>& ref=*this;
	reel res(0);
	for (int i=0;i<M;i++)
		for (int j=0;j<N;j++)
			res+=std::abs(ref(i,j));
	return res;
}
template <typename T,int SIZE,int M,int N> inline reel Matrix_Base<T,SIZE,M,N>::normL2()const
{
	return norm();
}
template <typename T,int SIZE,int M,int N> template <typename TT,int SSIZE> inline reel Matrix_Base<T,SIZE,M,N>::dot(const Matrix_Base<TT,SSIZE,M,N>& Two)const
{
	const Matrix_Base<T,SIZE,M,N>& ref=*this;
	reel res(0);
	for (int i=0;i<M;i++)
		for (int j=0;j<N;j++)
			res+=ref(i,j)*Two(i,j);
	return res;
}
template <typename T,int SIZE,int M,int N> template <typename TT,int SSIZE,int NN> inline typename Matrix_Product_Scheme_Trait< typename Matrix_Base<T,SIZE,M,N>::self, typename Matrix_Base<TT,SSIZE,N,NN>::self >::result_type Matrix_Base<T,SIZE,M,N>::operator*(const Matrix_Base<TT,SSIZE,N,NN>& M1) const {
	typename Matrix_Product_Scheme_Trait< typename Matrix_Base<T,SIZE,M,N>::self, typename Matrix_Base<TT,SSIZE,N,NN>::self >::result_type temp;
	const Matrix_Base<T,SIZE,M,N>&  ref=*this;
	return temp.prod(ref,M1);
}
template <typename T,int SIZE,int M,int N> template <typename TT,int SSIZE> inline typename Matrix_Base<T,SIZE,M,N>::self& Matrix_Base<T,SIZE,M,N>::mult(const Matrix_Base<TT,SSIZE,M,N>& Two) {
	typename Matrix_Strictly_Additive_Scheme_Trait< typename Matrix_Base<T,SIZE,M,N>::self, typename Matrix_Base<TT,SSIZE,M,N>::self >::result_type& ref=Child();
	return ref.mult(Two);
}
template <typename T,int SIZE,int M,int N> template <typename TT,int SSIZE> inline typename Matrix_Base<T,SIZE,M,N>::self& Matrix_Base<T,SIZE,M,N>::div(const Matrix_Base<TT,SSIZE,M,N>& Two) {
	typename Matrix_Strictly_Additive_Scheme_Trait< typename Matrix_Base<T,SIZE,M,N>::self, typename Matrix_Base<TT,SSIZE,M,N>::self >::result_type& ref=Child();
	return ref.div(Two);
}
template <typename T,int SIZE,int M,int N> template <typename TT,int SSIZE> inline bool Matrix_Base<T,SIZE,M,N>::operator==(const Matrix_Base<TT,SSIZE,M,N>& Two) const {
	const Matrix_Base<T,SIZE,M,N>& ref=*this;
	bool res=true;
	for (int i=0;(i<M)&&(res);i++)
		for (int j=0;(j<N)&&(res);j++)
			res=(ref(i,j)==Two(i,j));
	return res;
}
template <typename T,int SIZE,int M,int N> template <typename TT,int SSIZE> inline typename Matrix_Additive_Scheme_Trait< typename Matrix_Base<T,SIZE,M,N>::self, typename Matrix_Base<TT,SSIZE,M,N>::self >::result_type Matrix_Base<T,SIZE,M,N>::operator+(const Matrix_Base<TT,SSIZE,M,N>& Two) const {
	typename Matrix_Additive_Scheme_Trait< typename Matrix_Base<T,SIZE,M,N>::self, typename Matrix_Base<TT,SSIZE,M,N>::self >::result_type temp;
	const Matrix_Base<T,SIZE,M,N>&  ref=*this;
	return temp.add(ref,Two);
}
template <typename T,int SIZE,int M,int N> template <typename TT,int SSIZE> inline typename Matrix_Additive_Scheme_Trait< typename Matrix_Base<T,SIZE,M,N>::self, typename Matrix_Base<TT,SSIZE,M,N>::self >::result_type Matrix_Base<T,SIZE,M,N>::operator-(const Matrix_Base<TT,SSIZE,M,N>& Two) const {
	typename Matrix_Additive_Scheme_Trait< typename Matrix_Base<T,SIZE,M,N>::self, typename Matrix_Base<TT,SSIZE,M,N>::self >::result_type temp;
	const Matrix_Base<T,SIZE,M,N>&  ref=*this;
	return temp.subst(ref,Two);
}
template <typename T,int SIZE,int M,int N> template <typename TT,int SSIZE> inline typename Matrix_Base<T,SIZE,M,N>::self& Matrix_Base<T,SIZE,M,N>::operator+=(const Matrix_Base<TT,SSIZE,M,N>& Two) {
	typename Matrix_Strictly_Additive_Scheme_Trait< typename Matrix_Base<T,SIZE,M,N>::self, typename Matrix_Base<TT,SSIZE,M,N>::self >::result_type& ref=Child();
	return ref.add(Two);
}
template <typename T,int SIZE,int M,int N> template <typename TT,int SSIZE> inline typename Matrix_Base<T,SIZE,M,N>::self& Matrix_Base<T,SIZE,M,N>::operator-=(const Matrix_Base<TT,SSIZE,M,N>& Two) {
	typename Matrix_Strictly_Additive_Scheme_Trait< typename Matrix_Base<T,SIZE,M,N>::self, typename Matrix_Base<TT,SSIZE,M,N>::self >::result_type& ref=Child();
	return ref.subst(Two);
}
template <typename T,int SIZE,int M,int N> inline void Matrix_Base<T,SIZE,M,N>::write(std::ostream& f) const {
	f.write((const char*) Tab,SIZE*sizeof(reel));
}

template <typename T,int SIZE,int M,int N> inline void Matrix_Base<T,SIZE,M,N>::read(std::istream& f) {
	f.read((char*) Tab,SIZE*sizeof(reel));
}

template <typename T1,int SIZE1,typename T2,int SIZE2,int M,int N> inline typename Matrix_Additive_Scheme_Trait< typename Matrix_Base<T1,SIZE1,M,N>::self ,typename Matrix_Base<T2,SIZE2,M,N>::self >::result_type mult(const Matrix_Base<T1,SIZE1,M,N> &M1,const Matrix_Base<T2,SIZE2,M,N> &M2)
{
	typename Matrix_Additive_Scheme_Trait< typename Matrix_Base<T1,SIZE1,M,N>::self ,typename Matrix_Base<T2,SIZE2,M,N>::self >::result_type res(M1);
	res.mult(M2);
	return res;
}
template <typename T1,int SIZE1,typename T2,int SIZE2,int M,int N> inline typename Matrix_Additive_Scheme_Trait< typename Matrix_Base<T1,SIZE1,M,N>::self ,typename Matrix_Base<T2,SIZE2,M,N>::self >::result_type div(const Matrix_Base<T1,SIZE1,M,N> &M1,const Matrix_Base<T2,SIZE2,M,N> &M2) 
{
	typename Matrix_Additive_Scheme_Trait< typename Matrix_Base<T1,SIZE1,M,N>::self ,typename Matrix_Base<T2,SIZE2,M,N>::self >::result_type res(M1);
	res.div(M2);
	return res;
}
template <typename T,int SIZE,int M,int N> inline typename Matrix_Base<T,SIZE,M,N>::self sqrt(const Matrix_Base<T,SIZE,M,N> &M1) 
{
	typename Matrix_Base<T,SIZE,M,N>::self res;
	for (int i=0;i<SIZE;i++) {
		assert(M1[i]>=0);
		res[i]=::sqrt(M1[i]);
	}
	return res;
}
template <typename T,int SIZE,int M,int N> inline typename Matrix_Base<T,SIZE,M,N>::self operator+(const reel& value,const Matrix_Base<T,SIZE,M,N>& M1) {
	return M1+value;
}
template <typename T,int SIZE,int M,int N> inline typename Matrix_Base<T,SIZE,M,N>::self operator-(const reel& value,const Matrix_Base<T,SIZE,M,N>& M1) {
	return typename Matrix_Base<T,SIZE,M,N>::self(value)-M1;
}
template <typename T,int SIZE,int M,int N> inline typename Matrix_Base<T,SIZE,M,N>::self operator*(const reel& value,const Matrix_Base<T,SIZE,M,N>& M1) {
	return M1*value;
}

NAMESPACE_END

#endif
